
%% Do likelihood ratio test of full model vs. model with no price or quality

clear
cd ../china/data
load hs10_indlist.raw;
load trade_share.raw;
indN=length(hs10_indlist);

ZZ = 50
log_likelihood_full=zeros(ZZ,1);
log_likelihood_restricted=zeros(ZZ,1);

reg_bound=43;

sample_select=zeros(51,1);
for y=1:50
if (y>1 & y<5)
sample_select(y,1)=1;
else
if (y>5 & y<13)
sample_select(y,1)=1;
else 
if (y>13 & y<43)
sample_select(y,1)=1;
end
end
end
end

for w=1:ZZ

%if (w>1 & w<5) | (w>5 & w<reg_bound)
%if sample_select(w,1)==1
if (w>1 & w<5) | (w>5 & w<43)


folder=sprintf('../china/estimation/hs10/regular/raw_%d',hs10_indlist(w,1))

cd (folder)


load data2_q_regular.mat
load beta_newq_regular.mat;
load EV_newq.mat

else

folder=sprintf('../china/estimation/hs10/biggest/raw_%d',hs10_indlist(w,1))

cd (folder)

load data2_q_biggest.mat
load beta_biggest.mat
load EV.mat

end

nExp = X;
nImp = M;
S=N*nExp*nExp;
delta=0.975;

failnum=0;

U=beta_p_sol*ep + beta_x_sol*switched + beta_c_sol*city_switched + xi_tilde_sol*lambda;
CSVF= U + delta*EV_sol;

index=ones(nExp,1)';
for i=1:nExp
    index(1,i)=N*nExp*(i-1)+1;
end

index=repmat(index,N,1);
index=index(:);
index=repmat(index,nExp,1);

% index is a vector that tells us what the first number of EV calculation
% is for each element of EV.

%TransProbs=rand(N*N*nExp*nExp); Too big!

A=zeros(nExp,N,S);

for i=1:S
    B = (index(i,1):index(i,1)+N*nExp-1)';
    B = reshape(B,N,nExp);
    A(:,:,i) = B';
end
% Each A is a (nExp x N) matrix arranged so as to calculate EV
% For example EV(1) is a double sum: sum of  sum of CSVF(1),(6),(11);
% sum of CSVF(2),(7),(12); etc.

Q2 = zeros(nExp, S);
for j=1:nExp
    Q = A(:,:,N*(j-1)+1);
    Q2 (:, (j-1)*nExp*N +1: j*nExp*N ) = repmat(Q,1,nExp);
end
% Q2 is the matrix for adding up to calculate the denominator of P

P = exp(CSVF) ./ sum( exp( CSVF(Q2) ) )' ;
log_likelihood_full(w,1) = sum ( log ( P ( impchoicestate) ) );

end

%%%%%%%%%%%%%%%%%%%%%



for w=1:ZZ

%if (w>1 & w<5) | (w>5 & w<reg_bound)
%if sample_select(w,1)==1
if (w>1 & w<5) | (w>5 & w<43)


folder=sprintf('/projects/programs/monarch/china/china_updated/second_round/estimation/hs10/regular/raw_%d',hs10_indlist_finalparam(w,1))

cd (folder)


load data2_q_regular_nopq.mat
load beta_newq_regular_nopq.mat;
load EV_newq_regular_nopq.mat

else

folder=sprintf('/projects/programs/monarch/china/china_updated/second_round/estimation/hs10/biggest/raw_%d',hs10_indlist_finalparam(w,1))

cd (folder)

load data2_q_biggest_nopq.mat
load beta_newq_biggest_nopq.mat
load EV_newq_nopq.mat

end


N=1;
nExp = X;
nImp = M;
%S=N*nExp*nExp;
S=nExp*nExp;
delta=0.975;

failnum=0;


U=beta_x_sol*switched + beta_c_sol*city_switched ;
CSVF= U + delta*EV_sol;

index=ones(nExp,1)';
for i=1:nExp
    index(1,i)=N*nExp*(i-1)+1;
end

index=repmat(index,N,1);
index=index(:);
index=repmat(index,nExp,1);

% index is a vector that tells us what the first number of EV calculation
% is for each element of EV.

%TransProbs=rand(N*N*nExp*nExp); Too big!

A=zeros(nExp,N,S);

for i=1:S
    B = (index(i,1):index(i,1)+N*nExp-1)';
    B = reshape(B,N,nExp);
    A(:,:,i) = B';
end
% Each A is a (nExp x N) matrix arranged so as to calculate EV
% For example EV(1) is a double sum: sum of  sum of CSVF(1),(6),(11);
% sum of CSVF(2),(7),(12); etc.

Q2 = zeros(nExp, S);
for j=1:nExp
    Q = A(:,:,N*(j-1)+1);
    Q2 (:, (j-1)*nExp*N +1: j*nExp*N ) = repmat(Q,1,nExp);
end
% Q2 is the matrix for adding up to calculate the denominator of P

P = exp(CSVF) ./ sum( exp( CSVF(Q2) ) )' ;

log_likelihood_restricted(w,1) = sum ( log ( P ( impchoicestate) ) );

end

lr_stat=2*(log_likelihood_full-log_likelihood_restricted);
trade_share=trade_share(1:ZZ);
trade_share=trade_share/sum(trade_share);
weighted_lr_stat=sum(lr_stat.*trade_share)

lr_stat_g3=zeros(ZZ,1);
for i=1:ZZ
  if lr_stat(i,1)>3
    lr_stat_g3(i,1)=1;
  end 
end
